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PREFACE 


In  December  1988,  the  Norfolk  District  Corps  of  Engineers  requested  that 
the  US  Army  Engineer  Waterways  Experiment  Station  (WES)  conduct  an  investiga¬ 
tion  to  assess  general  changes  in  circulation,  currents,  and  sedimentation 
associated  with  three  proposed  overdeepening  plans  for  the  Newport  News  Chan¬ 
nel  between  the  1-664  Bridge-Tunnel  crossing  and  the  Hampton  Roads 
Bridge-Tunnel  crossing. 

The  study  was  conducted  by  personnel  of  the  Hydraulics  Laboratory,  WES, 
under  the  general  direction  of  Messrs.  F.  A.  Herrmann,  Jr.,  Chief  of  the  Hy¬ 
draulics  Laboratory;  R.  A.  Sager,  Assistant  Chief  of  the  Hydraulics  Labora¬ 
tory;  W.  H.  McAnally,  Jr.,  Chief  of  the  Estuaries  Division;  and  W.  D.  Martin, 
Chief  of  the  Estuarine  Engineering  Branch.  The  study  was  conducted  by  Dr. 
Hsin-Chi  J.  Lin,  with  technical  consultation  supplied  by  Messrs.  S.  B.  Heltzel 
and  M.  A.  Granat,  all  of  the  Estuarine  Engineering  Branch.  This  report  was 
prepared  by  Mr.  Martin  and  Dr.  Lin  and  edited  by  Mrs.  Marsha  C.  Gay  of  the 
Information  Technology  Laboratory,  WES. 

Acting  Commander  and  Director  of  WES  during  preparation  of  this  report 
was  LTC  Jack  R.  Stephens,  EN.  Technical  Director  was  Dr.  Robert  W.  Whalln. 
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CONVERSION  FACTORS,  NON-SI  TO  SI  (METRIC) 
UNITS  OF  MEASUREMENT 


Non-SI  units  of  measurement 
(metric)  units  as  follows: 

Multiply 
cubic  feet 
feet 

pounds  (force) - 
second  per 
foot  per  foot 


used  In  this  report  can 

By 

0.02831685 

0.3048 

47.88026 


be  converted  to  SI 

To  Obtain 
cubic  metres 
metres 

pascals-second 
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NEWPORT  NEWS  CHANNEL  DEEPENING  STUDY,  VIRGINIA 


Numerical  Model  Investigation 


PART  I:  INTRODUCTION 

Background 

1.  The  deepening  of  the  Newport  News  Channel  from  45  to  55  ft*  was 
studied  by  the  US  Army  Engineer  Waterways  Experiment  Station  (WES)  as  a  por¬ 
tion  of  the  Norfolk  Harbor  and  Channels  deepening  study  (Richards  and  Morton 
1983).  That  study  was  conducted  on  the  Chesapeake  Bay  physical  model  located 
on  Kent  Island  in  Matapeake,  MD. 

2.  The  effects  on  sedimentation  of  this  deepening  were  also  investi¬ 
gated  by  WES  using  the  hybrid  modeling  approach.  This  approach  combined  the 
physical  model  results  with  a  numerical  analysis  using  the  WES  TABS-2  system 
of  numerical  models  (Berger  et  al.  1985). 

3.  Other  WES  studies  in  this  area  evaluated  the  effects  of  the  1-664 
Bridge-Tunnel  crossing  (Heltzel  1988)  and  the  enlargement  of  the  Craney  Island 
disposal  area  (Heltzel  and  Granat  1988  and  Bottin  1984) . 

4.  To  expand  the  Craney  Island  facility,  a  levee  is  to  be  constructed 
using  expansion  Plan  B  (Heltzel  and  Granat  1988) .  The  Norfolk  District  Corps 
of  Engineers  wished  to  investigate  obtaining  the  material  to  construct  the 
levee  by  overdeepening  the  Newport  News  Channel  between  the  1-664  crossing  and 
the  Hampton  Roads  Bridge-Tunnel  crossing  (Figure  1)  from  the  current  depth  of 
55  ft.  Three  alternative  depths  of  57,  64,  and  70  ft  were  evaluated. 

Purpose 


5.  The  objective  of  this  study  was  to  use  available  numerical  models  to 
evaluate  general  changes  in  circulation,  currents,  and  sedimentation  associ¬ 
ated  with  the  overdeepening  of  the  Newport  News  Channel.  Additionally,  the 
effects  of  the  overdeepening  on  the  reported  estuarine  circulation  cell  (flow 


*  A  table  of  factors  for  converting  non-SI  units  of  measurement  to  SI 
(metric)  units  of  measurement  is  found  on  page  3. 
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convergence)  o££  Hampton  Flats  and  Newport  News  Point  was  to  be  assessed. 


Scope 


6.  The  numerical  modeling  was  designed  to  evaluate  relative  changes  In 
hydrodynamics  and  sedimentation  adjacent  to  the  study  area.  The  sedimentation 
comparisons  £ocused  on  two  areas  adjacent  to  the  Newport  News  Channel  and  the 
deepened  portion  o£  the  channel.  The  c££-channel  areas  were  designated  A  and 
B  and  were  the  same  as  those  reported  by  Heltzel  and  Granat  (1988).  These 
zones  are  shown  In  Figure  2.  In  these  relatively  low  velocity  areas,  the 
sediment  study  £ocused  on  changes  In  cohesive  sediment  transport.  The  channel 
comparisons  were  made  In  an  area  o£  relatively  high  velocities  and  £ocused  on 
noncoheslve  sediment  transport. 

7.  The  circulation  cell  o££  Newport  News  Point  and  Hampton  Flats  was 
addressed  by  comparing  study  results  with  previously  compiled  data  on  the 
phenomenon  (Heltzel  and  Granat  1988) . 
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Figure  2.  Critical  sedimentation  zones  of  interest 


PART  II:  NUMERICAL  MODELING  APPROACH 


The  Ntunerlcal  Models 


8.  The  Corps  numerical  modeling  system,  Open-Channel  Flow  and  Sedimen¬ 
tation,  TABS-2  (Thomas  and  McAnally  1985),  was  used  in  this  investigation. 

The  two  primary  finite  element  numerical  codes  used  were  A  Two-Dimensional 
Model  for  Free  Surface  Flows  (RMA-2V)  and  Sediment  Transport  in  Unsteady  Two- 
Dimensional  Flows,  Horizontal  Plane  (STUDH).  Both  computer  codes  employ  the 
finite  element  method  to  solve  the  depth-integrated  governing  equations.  A 
description  of  RMA-2V  and  STUDH  appears  in  Appendix  A. 

Newport  News  Channel  Computational  Meshes 

9.  The  computational  mesh  used  in  this  study  was  a  modified  version  of 
the  mesh  used  in  the  Lower  James  River  circulation  study  (Heltzel  and  Granat 
1988).  The  following  modifications  were  included: 

a.  The  manner  in  which  the  I-66A  North  and  South  Islands  were 
represented  in  the  model  was  revised.  They  were  previously 
represented  by  elements  with  Increased  roughness  coefficients. 
For  this  study,  they  were  modeled  as  solid  structures  with  slip 
flow  boundaries. 

b.  The  new  small  boat  harbor  at  Newport  News  Point  was  added  to  the 
mesh. 

c.  The  mesh  resolution  was  Increased  in  the  vicinity  of  the  North 
and  South  islands  and  in  the  Hampton  Flats  area. 

d.  Additional  resolution  was  added  to  allow  for  the  modeling  of  the 
overdeepened  Newport  News  Channel. 

The  limits  of  the  overdeepened  reach  are  shown  in  Figure  1  .  An  easement  ex¬ 
tending  500  ft  on  either  side  of  the  1-664  tunnel  crossing  was  not  over¬ 
deepened. 

10.  The  revised  mesh,  shown  in  Figure  3,  contains  2,672  nodes  and  933 
elements.  This  mesh  Incorporated  as  base  conditions  those  tested  as  expansion 
Plan  B  for  Craney  Island  in  the  Lower  James  River  circulation  study  (Heltzel 
and  Granat  1988)  with  the  exception  of  the  addition  of  the  small  boat  harbor 
adjacent  to  the  Island  north  of  the  1-664  crossing.  A  detail  of  the  Craney 
Island  expansion  configuration  incorporated  in  the  mesh  is  shown  in  Figure  4. 
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Figure  4.  Craney  Island  expansion  plan  B  used  In  Newport  News  Channel 

deepening  study 


Figure  5  shows  a  detail  of  the  computational  mesh  in  the  vicinity  of  the  1-664 
bridge  crossing  and  Hampton  Flats. 


Testing  Conditions 


11.  Three  plans  for  overdeepening  the  channel  from  55  ft  (Plan  0)  were 
tested,  representing  depths  of  57  ft  (Plan  1),  64  ft  (Plan  2),  and  70  ft 
(Plan  3). 

12.  Boundary  conditions  were  Identical  to  those  used  in  the  Lower  James 
River  circulation  study  (Heltzel  and  Granat  1988).  These  were  developed  from 
physical  model  data  collected  In  the  Chesapeake  Bay  physical  model  and  used  In 
the  Norfolk  Harbor  and  Channels  deepening  study  (Richards  and  Morton  1983). 

The  lower  or  bay  boundary  of  the  model  was  represented  by  water-surface  eleva¬ 
tion  data.  The  mean  range  tide  (2.5  ft  at  Old  Point  Comfort)  was  used. 
Depth-averaged  velocity  data  represented  the  Inflows  of  the  Elizabeth  and 
Upper  James  rivers.  The  long-term  average  James  River  discharge  of  8,900  cfs 
and  Elizabeth  River  discharge  of  300  cfs  were  used. 

13.  In  general,  the  hydrodynamic  and  sediment  coefficients  and  modeling 
procedures  used  in  the  Lower  James  River  circulation  study  were  used  In  this 
study  with  the  exception  of  the  eddy  viscosity  coefficients,  which  were  re¬ 
laxed  to  a  value  of  50  Ib-sec/ft/ft  in  the  vicinity  of  Newport  News  Point  to 
reproduce  the  circulation  cells  Indicated  by  physical  model  results.  The 
Plan  B  expansion  scheme  from  the  Lower  James  study  was  reactivated  and  tested 
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to  ensure  that  the  model  was  yielding  consistent  results.  Once  the  mesh  modi¬ 
fications  were  incorporated,  a  similar  check  indicated  a  slight  change  in  the 
magnitude  of  the  Plan  B  results  from  those  previously  reported  (Heltzel  and 
Granat  1988).  Maximum  differences  were  0.1  fps  or  less  in  areas  A  and  B  and 
0.3  fps  in  the  main  channel.  Therefore,  the  results  reported  herein  are  not 
directly  comparable  to  the  previous  study  for  small  changes. 
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PART  III;  MODELING  RESULTS  AND  DISCUSSION 


Hydrodynamic  Results 

14.  The  channel  overdeepening  resulted  in  no  discernible  change  in  the 
two-dimensional  circulation  patterns  in  the  Lower  James  River.  The  addition 
of  the  small  boat  harbor  and  the  increased  mesh  resolution  resulted  in  a  more 
clearly  defined  eddy  over  the  Hampton  Flats  off  Newport  News  Point  than  that 
previously  reported  (Heltzel  and  Granat  1988)  during  the  period  around  slack 
before  flood  (hours  15,  16,  and  17).  Detailed  vector  plots  for  hours  15,  16, 
17,  and  18  for  each  plan  are  shown  in  Plates  1-16. 

15.  Figures  3  and  5  show  the  location  of  seven  nodes  that  were  used  for 
comparison  of  the  three  deepening  plans.  Differences  in  the  magnitudes  of  the 
maximum  ebb  and  flood  velocities  for  the  base  and  the  three  plans  were  com¬ 
pared,  and  the  results  are  suua&arlzed  in  Table  1.  Node  18  is  located  in  the 
center  of  the  Willoughby  Bay  area  of  interest.  Nodes  90,  146,  and  198  are 
located  in  the  Hampton  Flats  area.  Nodes  1972,  1980,  and  2000  are  located  in 
the  Newport  News  Channel  proper. 

16.  As  can  be  seen  in  Table  1,  there  was  no  difference  in  the  velocity 
magnitudes  for  any  of  the  plans  tested  in  the  Willoughby  Bay  area.  Plans  1 
and  2  had  virtually  no  effect  on  the  velocities  in  the  Hampton  Flats  area. 

Flan  3  showed  a  slight  but  measurable  decrease  in  velocities  of  the  flats,  the 
maximum  decrease  being  less  than  0.1  fps  at  node  198.  Velocities  in  the  New¬ 
port  News  channel  were  uniformly  reduced,  as  would  be  expected,  by  the  channel 
overdeepening.  The  maximum  reduction  in  ebb  velocity  was  0.2  fps  at 

node  1972.  The  maximum  decrease  in  flood  velocity  was  0.1  fps  at  node  1972. 

17.  Plates  17  and  18  illustrate  the  time-history  plots  of  the  ebb  and 
flood  velocities  for  the  selected  nodes  comparing  Plan  0  (55-ft  depth)  with 
Plan  1  (57-ft  depth).  It  can  be  seen  that  these  plots  are  virtually  identi¬ 
cal.  Plates  19  and  20  illustrate  the  same  comparison  for  Plan  0  and  Plan  2 
(64-f t  depth) .  Only  slight  reductions  in  the  ebb  and  flood  velocities  were 
noted  for  nodes  1972  and  2000,  both  of  which  are  located  in  the  channel 
proper.  Plates  21  and  22  illustrate  the  comparison  for  Plan  0  and  Plan  3 
(70-ft  depth).  Again,  little  discernible  change  in  velocities  is  observed 
outside  the  channel  area,  while  slight  reductions  were  observed  at  nodes  1972, 
1980,  and  2000  in  the  channel. 
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18.  Water-surface  elevations  were  also  compared  for  the  base  and  three 
plan  deepenlngs.  These  tlme-hlstorles  were  Identical  for  each  comparison. 
Indicating  no  effect  on  water-surface  elevations  or  phasing  of  the  tide. 

These  plots  are,  therefore,  not  Included  In  this  report. 

Sediment  Results 

19.  The  approach  to  analyzing  the  sediment  results  consisted  of  quali¬ 
tative  comparisons  between  the  base  (Plan  0)  and  the  three  deepening  plans. 

The  procedure  followed  and  parameters  selected  duplicated  those  used  In  the 
1-664  Bridge-Tunnel  study  (Heltzel  1988). 

20.  The  areas  with  relatively  low  velocity  that  were  of  Interest  In 
Willoughby  Bay  and  Hampton  Flats  were  analyzed  using  cohesive  sediment  model¬ 
ing  techniques.  The  higher  velocity  areas  In  the  channel  proper  were  analyzed 
using  noncoheslve  sediment  modeling  techniques. 

21.  The  sediment  analysis  was  limited  to  three  areas,  area  A  (see  Fig¬ 
ure  2)  In  Willoughby  Bay,  area  B  In  Hampton  Flats,  and  the  overdeepened  por¬ 
tion  of  the  Newport  News  Channel.  For  reporting  purposes,  the  predicted 
shoaling  volume  for  each  element  In  the  zone  of  Interest  was  combined  to  pro¬ 
duce  a  cumulative  rate  for  each  area.  The  plan  rate  was  then  divided  by  the 
base  (Plan  0)  rate  to  produce  a  shoaling  index,  which  could  then  be  used  as  a 
basis  for  comparison  between  the  plans.  The  results  of  this  comparison  are 
summarized  In  the  following  tabulation.  Plans  1  and  2  showed  a  slight  in¬ 
crease  In  shoaling  in  the  Willoughby  Bay  area.  Plan  3  Indicated  a  6  percent 
Increase  In  shoaling.  All  plans  showed  a  decrease  In  shoaling  in  the  Hampton 
Flats  area,  with  the  maximum  decrease  being  7  percent.  The  Newport  News  Chan¬ 
nel  showed  very  slight  Increases  In  the  shoaling  rate,  the  Increases  amounting 
to  less  than  1  percent. 


Shoaling  Index* 


Area 

Plan  1 

Plan  2 

Plan  3 

A 

1.01 

1.03 

1.06 

B 

0,93 

0.94 

0.95 

Newport  News 

1.00+ 

1.00+ 

1.00+ 

Channel 


*  1.00+  Indicates  amount  less  than  1  percent. 
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22.  The  deepening  will  cause  a  redistribution  of  concentrations  of 
cohesive  sediments  that  will  account  for  generally  decreased  shoaling  In  the 
Hampton  Flats  area  and  Increased  shoaling  In  the  Willoughby  Bay  area.  It 
should  be  pointed  out  that  the  Increases  are  reported  as  percentages  and  rep¬ 
resent  relatively  small  absolute  Increases  In  shoaling.  Reduced  velocities  In 
the  channel  proper  and  a  small  but  definite  tendency  to  shoal  will  result  In  a 
gradual  filling  of  the  overdeepened  areas  and  an  eventual  return  to  base 
conditions. 
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PART  IV:  ESTUARINE  CIRCULATION  AND  FLOW  CONVERGENCE: 
HAMPTON  FLATS  AND  NEWPORT  NEWS  POINT 


23.  A  previous  WES  report  (Heltzel  and  Granat  1988)  summarized  the  In¬ 
formation  available  regarding  the  estuarine  circulation  and  flow  convergence 
observed  off  Newport  News  Point  and  the  Hampton  Flats.  This  summary  Included 
Information  from  previous  model  studies,  both  physical  and  numerical  (Brogdon 
and  Bobb  1967,*  Heltzel  1988),  data  collected  by  the  Virginia  Institute  for 
Marine  Sciences  (VIMS)  (Byrne  et  al.  1987),  and  WES  field  data  collected  In 
1986. 

24.  The  potential  Impacts  to  the  reported  frontal  system,  a  three- 
dimensional  phenomenon,  cannot  be  quantified  with  Information  from  either  past 
or  present  studies.  However,  Inferences  can  be  drawn  from  the  present 
two-dimensional  study. 

25.  Based  on  the  hydrodynamic  results,  there  will  be  no  discernible 
changes  In  two-dimensional  circulation  patterns  In  the  area  off  Newport  News 
Point  or  over  the  Hampton  Flats  due  to  the  channel  overdeepening.  Velocity 
magnitudes  will  decrease  by  such  a  small  amount  that  the  Impacts  of  this 
change  should  also  be  negligible.  The  tide  phasing  and  elevations  were  also 
unaffected  by  the  deepening. 


*  N.  J.  Brogdon,  Jr.,  and  W.  H.  Bobb.  1967.  "Effects  of  Proposed  Waterfront 
Developments  at  Newport  News  Point  on  Tides,  Currents,  Salinities,  and 
Shoaling,"  Draft  Report,  US  Army  Engineer  Waterways  Experiment  Station, 
Vicksburg,  MS. 
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PART  V:  CONCLUSIONS 


26.  Comparisons  of  base  and  channel  velocities  for  overdeepening  plans 
Indicate  no  velocity  Increases  In  the  Willoughby  Bay  area.  The  maximum  veloc¬ 
ity  change  Identified  In  the  Hampton  Flats  area  was  a  less  than  0.1-fps 
decrease  In  the  maximum  flood  tide  velocity.  Velocities  In  the  channel  proper 
decreased  a  maximum  of  0.2  fps. 

27.  No  circulation  changes  were  Identified  In  base-to-plan  comparisons 
of  vector  plots.  Additionally,  no  change  In  tidal  phasing  or  water-surface 
elevations  was  detected. 

28.  Qualitatively,  the  various  deepening  plans  resulted  In  a  redistri¬ 
bution  of  cohesive  sediments  with  a  net  loss  over  the  Hampton  Flats  area  and  a 
net  Increase  In  the  Willoughby  Bay  area.  The  channel  will  experience  a  slow 
rate  of  shoaling  due  to  the  overdeepening.  However,  all  changes  In  sedimen¬ 
tation  were  small  In  absolute  volumes. 

29.  The  formation  and  location  of  the  two-dimensional  circulation  cell 
off  Newport  News  Point  were  unaffected  by  any  of  the  plans  addressed. 

30.  The  reported  frontal  effect  off  Newport  News  Point  Is  a  three- 
dimensional  density  current-driven  phenomenon  and  cannot  be  quantified  within 
the  scope  of  this  two-dimensional  analysis.  However,  no  evidence  was  gener¬ 
ated  by  this  study  that  would  Indicate  that  the  channel  overdeepening  will 
affect  the  frontal  formation  or  propagation. 
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Table  1 

Maximum  Velocity  Changes  at  Selected  Nodes,  fps 
(Plan  Minus  Base) 


Area 

Node  No. 

Plan  0 
Velocity 

Chang 
Plan  1 

e ,  Plan  Minus 
Plan  2 

Base 

Plan  3 

Ebb 

Flood 

Ebb 

Flood 

Ebb 

Flood 

Ebb 

Flood 

Willoughby 

18 

0.77 

0.75 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

Bay  (area  A) 

Hampton  Flats 

90 

0.58 

0.62 

0.0 

0.0 

0.0 

T+ 

0.0 

T+ 

(area  B) 

146 

0.64 

0.60 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

198 

0.83 

0.85 

0.0 

0.0 

T- 

T- 

T- 

T- 

Newport  News 

1972 

2.53 

1.94 

-0.1 

T- 

-0.1 

-0.1 

-0.2 

-0.1 

Channel 

1980 

1.80 

1.58 

0.0 

0.0 

T- 

T- 

T- 

T- 

2000 

1.20 

1.07 

0.0 

T- 

T- 

T- 

T- 

-0.1 

Note:  Velocities  rounded  to  nearest  0.1  fps. 

T+  or  T-  indicates  values  less  than  0.1  fps 


PLATE  2 


VELOCITY  VECTOR 


DEEPENING  PLAN  2 
(64-FT  DEPTH) 
EXPANSION  PLAN  B 
HOUR  15 


PLATE  9 


PLATE  17 


VELOCITY  TIME-HISTORIES 
PLAN  0  VERSUS  PLAN  1 
NODES  1972-2000 


PLATE  19 


PLATE  21 


VELOCITY  TIME-HISTORIES 
PLAN  0  VERSUS  PLAN  3 
NODES  1972-2000 


APPENDIX  A:  THE  TABS-2  SYSTEM 


1.  TABS-2  Is  a  collection  of  generalized  computer  programs  and  utility 
codes  integrated  into  a  numerical  modeling  system  for  studying  two-dimensional 
hydrodynamics,  sedimentation,  and  transport  problems  in  rivers,  reservoirs, 
bays,  and  estuaries.  A  schematic  representation  of  the  system  is  shown  in 
Figure  A1 .  It  can  be  used  either  as  a  stand-alone  solution  technique  or  as  a 
step  in  the  hybrid  modeling  approach.  The  basic  concept  is  to  calculate 
water-surface  elevations,  current  patterns,  sediment  erosion,  transport  and 
deposition,  the  resulting  bed  surface  elevations,  and  the  feedback  to  hydrau¬ 
lics.  Existing  and  proposed  geometry  can  be  analyzed  to  determine  the  Impact 
on  sedimentation  of  project  designs  and  to  determine  the  impact  of  project 
designs  on  salinity  and  on  the  stream  system.  The  system  is  described  in  de¬ 
tail  by  Thomas  and  McAnally  (1985). 

2.  The  three  basic  components  of  the  system  are  as  follows: 

a.  "A  Two-Dimensional  Model  for  Free  Surface  Flows,"  RMA-2V. 

b.  "Sediment  Transport  in  Unsteady  2-Dimensional  Flows,  Horizontal 
Plane,"  STUDH. 

£.  "Two-Dimensional  Finite  Element  Program  for  Water  Quality," 
RMA-4. 

3.  RMA-2V  is  a  finite  element  solution  of  the  Reynolds  form  of  the 
Navier-Stokes  equations  for  turbulent  flows.  Friction  is  calculated  with 
Manning’s  equation  and  eddy  viscosity  coefficients  are  used  to  define  the 
turbulent  losses.  A  velocity  form  of  the  basic  equation  is  used  with  side 
boundaries  treated  as  either  slip  or  static.  The  model  automatically  recog¬ 
nizes  dry  elements  and  corrects  the  mesh  accordingly.  Boundary  conditions  may 
be  water-surface  elevations,  velocities,  or  discharges  and  may  occur  inside 
the  mesh  as  well  as  along  the  edges. 

TWS-2 


Figure  Al.  TABS-2  schematic 
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4.  The  sedimentation  model,  STUDH,  solves  the  convection-diffusion 
equation  with  bed  source  terms.  These  terms  are  structured  for  either  sand  or 
cohesive  sediments.  The  Ackers-White  (1973)  procedure  is  used  to  calculate  a 
sediment  transport  potential  for  the  sands  from  which  the  actual  transport  is 
calculated  based  on  availability.  Clay  erosion  is  based  on  work  by  Parthen- 
iades  (1962)  and  Ariathurai  and  the  deposition  of  clay  utilizes  Krone's  equa¬ 
tions  (Ariathural,  MacArthur,  and  Krone  1977).  Deposited  material  forms 
layers,  as  shown  in  Figure  A2,  and  bookkeeping  allows  up  to  10  layers  at  each 
node  for  maintaining  separate  material  types,  deposit  thickness,  and  age.  The 
code  uses  the  same  mesh  as  RMA-2V. 

5.  Salinity  calculations,  RMA-4,  are  made  with  a  form  of  the 
convective-diffusion  equation  which  has  general  source-sink  terms.  Up  to 
seven  conservative  substances  or  substances  requiring  a  decay  term  can  be 
routed.  The  code  uses  the  same  mesh  as  RMA-2V. 

6.  Each  of  these  generalized  computer  codes  can  be  used  as  a  stand¬ 
alone  program,  but  to  facilitate  the  preparation  of  input  data  and  to  aid  in 
analyzing  results,  a  family  of  utility  programs  was  developed  for  the  follow¬ 
ing  purposes: 

a.  Digitizing 

b.  Mesh  generation 

c.  Spatial  data  management 

d.  Graphical  output 

e.  Output  analysis 

£.  File  management 

£.  Interfaces 

h.  Job  control  language 

Finite  Element  Modeling 


7.  The  TABS-2  numerical  models  used  in  this  effort  employ  the  finite 
element  method  to  solve  the  governing  equations.  To  help  those  who  are  un¬ 
familiar  with  the  method  to  better  understand  this  report,  a  brief  description 
of  the  method  is  given  here. 

8.  The  finite  element  method  approximates  a  solution  to  equations  by 
dividing  the  area  of  interest  into  smaller  subareas,  which  are  called  ele¬ 
ments.  The  dependent  variables  (e.g.,  water-surface  elevations  and  sediment 
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concentrations)  are  approximated  over  each  element  by  continuous  functions 
which  interpolate  in  terms  of  unknown  point  (node)  values  of  the  variables. 

An  error,  defined  as  the  deviation  of  the  approximation  solution  from  the  cor¬ 
rect  solution,  is  minimized.  Then,  when  boundary  conditions  are  imposed,  a 
set  of  solvable  simultaneous  equations  is  created.  The  solution  is  continuous 
over  the  area  of  interest. 

9.  In  one-dimensional  problems,  elements  are  line  segments.  In  two- 
dimensional  problems,  the  elements  are  polygons,  usually  either  triangles  or 
quadrilaterals.  Nodes  are  located  on  the  edges  of  elements  and  occasionally 
inside  the  elements.  The  interpolating  functions  may  be  linear  or  higher 
order  polynomials.  Figure  A2  Illustrates  a  quadrilateral  element  with  eight 
nodes  and  a  linear  solution  surface  where  F  is  the  interpolating  function. 

10.  Most  water  resource  applications  of  the  finite  element  method  use 
the  Galerkln  method  of  weighted  residuals  to  minimize  error.  In  this  method 
the  residual,  the  total  error  between  the  approximate  and  correct  solutions, 
is  weighted  by  a  function  that  is  identical  with  the  interpolating  function 
and  then  minimized.  Minimization  results  in  a  set  of  simultaneous  equations 
in  terms  of  nodal  values  of  the  dependent  variable  (e.g.  water-surface  eleva¬ 
tions  or  sediment  concentration).  The  time  portion  of  time-dependent  problems 
can  be  solved  by  the  finite  element  method,  but  it  is  generally  more  efficient 
to  express  derivatives  with  respect  to  time  in  finite  difference  form. 

The  Hydrodynamic  Model,  RMA-2V 


Applications 

11.  This  program  is  designed  for  far-field  problems  in  which  vertical 
accelerations  are  negligible  and  the  velocity  vectors  at  a  node  generally 
point  in  the  same  directions  over  the  entire  depth  of  the  water  column  at  any 
instant  of  time.  It  expects  a  homogeneous  fluid  with  a  free  surface.  Both 
steady  and  unsteady  state  problems  can  be  analyzed.  A  surface  wind  stress  can 
be  imposed. 

12.  The  program  has  been  applied  to  calculate  flow  distribution  around 
Islands;  flow  at  bridges  having  one  or  more  relief  openings,  in  contracting 
and  expanding  reaches,  into  and  out  of  off-channel  hydropower  plants,  at  river 
junctions,  and  into  and  out  of  pumping  plant  channels;  and  general  flow  pat¬ 
terns  in  rivers,  reservoirs,  and  estuaries. 
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Limitations 


13.  This  program  Is  not  designed  for  near-fleld  problems  where  flow- 
structure  interactions  (such  as  vortices,  vibrations,  or  vertical  accelera¬ 
tions)  are  of  interest.  Areas  of  vertically  stratified  flow  are  beyond  this 
program's  capability  unless  it  is  used  in  a  hybrid  modeling  approach.  It  is 
two-dimensional  in  the  horizontal  plane,  and  zones  where  the  bottom  current  is 
in  a  different  direction  from  the  surface  current  must  be  analyzed  with  con¬ 
siderable  subjective  judgement  regarding  long-term  energy  considerations.  It 
is  a  free-surface  calculation  for  subcritical  flow  problems. 


Governing  equations 

14.  The  generalized  computer  program  RMA-2V  solves  the  depth-integrated 


equations  of  fluid  mass  and  momentum  conservation  in  two  horizontal  direc¬ 


tions.  The  form  of  the  solved  equations  is 


where 


h  =  depth 

u,v  =  velocities  in  the  Cartesian  directions 
x,y,t  =  Cartesian  coordinates  and  time 
p  =  density 


(Al) 


(A2) 


(A3) 
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e  =  eddy  viscosity  coefficient,  for  xx  =  normal  direction  on 

x-axis  surface;  yy  =  normal  direction  on  y-axis  surface;  xy 
and  yx  =  shear  direction  on  each  surface 

g  =  acceleration  due  to  gravity 

a  =  elevation  of  bottom 

n  =  Manning's  n  value 

1.486  =  conversion  from  SI  (metric)  to  non-SI  units 
4  =  empirical  wind  shear  coefficient 
=  wind  speed 
tp  =  wind  direction 

w  =  rate  of  earth's  angular  rotation 
=  local  latitude 

15.  Equations  Al,  A2,  and  A3  are  solved  by  the  finite  element  method 
using  Galerkln  weighted  residuals.  The  elements  may  be  either  quadrilaterals 
or  triangles  and  may  have  curved  (parabolic)  sides.  The  shape  functions  are 
quadratic  for  flow  and  linear  for  depth.  Integration  in  space  is  performed  by 
Gaussian  integration.  Derivatives  in  time  are  replaced  by  a  nonlinear  finite 
difference  approximation.  Variables  are  assumed  to  vary  over  each  time  inter¬ 
val  in  the  form 


f(t)  =  ffO)  +  at  +  bt*^  tp  <  t  <  t  (A4) 

which  is  differentiated  with  respect  to  time,  and  cast  in  finite  difference 
form.  Letters  a  ,  b  ,  and  c  are  constants.  It  has  been  found  by  experi¬ 
ment  that  the  best  value  for  c  is  1.5  (Norton  and  King  1977). 

16.  The  solution  is  fully  implicit  and  the  set  of  simultaneous  equations 
is  solved  by  Newton-Raphson  iteration.  The  computer  code  executes  the  solu¬ 
tion  by  means  of  a  front-type  solver  that  assembles  a  portion  of  the  matrix 
and  solves  it  before  assembling  the  next  portion  of  the  matrix.  The  front 
solver's  efficiency  is  largely  independent  of  bandwidth  and  thus  does  not 
require  as  much  care  in  formation  of  the  computational  mesh  as  do  traditional 
solvers. 

17.  The  code  RMA-2V  is  based  on  the  earlier  version  RMA-2  (Norton  and 
King  1977)  but  differs  from  it  in  several  ways.  It  is  formulated  in  terms  of 
velocity  (v)  instead  of  unit  discharge  (vh) ,  which  improves  some  aspects  of 
the  code's  behavior;  it  permits  drying  and  wetting  of  areas  within  the  grid; 
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and  it  permits  specification  of  turbulent  exchange  coefficients  in  directions 
other  than  along  the  x-  and  z-axes.  For  a  more  complete  description,  see 
Appendix  F  of  Thomas  and  McAnally  (1985). 

The  Sediment  Transport  Model,  STUDH 


Applications 

18.  STUDH  can  be  applied  to  clay  and/or  sand  bed  sediments  where  flow 
velocities  can  be  considered  two-dimensional  (i.e.,  the  speed  and  direction 
can  be  satisfactorily  represented  as  a  depth-averaged  velocity) .  It  is  useful 
for  both  deposition  and  erosion  studies  and,  to  a  limited  extent,  for  stream 
width  studies.  The  program  treats  two  categories  of  sediment:  noncohesive, 
which  is  referred  to  as  sand  here,  and  cohesive,  which  is  referred  to  as  clay. 
Limitations 

19.  Both  clay  and  sand  may  be  analyzed,  but  the  model  considers  a 
single,  effective  grain  size  for  each  and  treats  each  separately.  Fall  veloc¬ 
ity  must  be  prescribed  along  with  the  water-surface  elevations,  x-velocity, 
y-velocity,  diffusion  coefficients,  bed  density,  critical  shear  stresses  for 
erosion,  erosion  rate  constants,  and  critical  shear  stress  for  deposition. 

20.  Many  applications  cannot  use  long  simulation  periods  because  of 
their  computation  cost.  Study  areas  should  be  made  as  small  as  possible  to 
avoid  an  excessive  number  of  elements  when  dynamic  runs  are  contemplated  yet 
must  be  large  enough  to  permit  proper  posing  of  boundary  conditions.  The  same 
computation  time  Interval  must  be  satisfactory  for  both  the  transverse  and 
longitudinal  flow  directions. 

21.  The  program  does  not  compute  water-surface  elevations  or  velocities; 
therefore  these  data  must  be  provided.  For  complicated  geometries,  the  numer¬ 
ical  model  for  hydrodynamic  computations,  RMA-2V,  is  used. 

Governing  equations 

22.  The  generalized  computer  program  STUDH  solves  the  depth-integrated 
convection-dispersion  equation  in  two  horizontal  dimensions  for  a  single  sedi¬ 
ment  constituent.  For  a  more  complete  description,  see  Appendix  G  of  Thomas 
and  McAnally  (1985).  The  form  of  the  solved  equation  is 


3C  ^  3C  ^  9C  3 

7—  +  UT-  +  V  —  =  -— 

3t  3x  3y  3x 


=  0 


(A5) 
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where 


C  =  concentration  of  sediment 

u  =  depth-integrated  velocity  in  x-direction 

V  =  depth-integrated  velocity  in  y-direction 
=  dispersion  coefficient  in  x-direction 
D  =  dispersion  coefficient  in  y-direction 

=  coefficient  of  concentration-dependent  source/sink  term 
=  coefficient  of  source/sink  term 

23.  The  source/sink  terms  in  Equation  A5  are  computed  in  routines  that 
treat  the  interaction  of  the  flow  and  the  bed.  Separate  sections  of  the  code 
handle  computations  for  clay  bed  and  sand  bed  problems. 

Sand  transport 

24.  The  source/sink  terms  are  evaluated  by  first  computing  a  potential 
sand  transport  capacity  for  the  specified  flow  conditions,  comparing  that 
capacity  with  the  amount  of  sand  actually  being  transported,  and  then  eroding 
from  or  depositing  to  the  bed  at  a  rate  that  would  approach  the  equilibrium 
value  after  sufficient  elapsed  time. 

25.  The  potential  sand  transport  capacity  in  the  model  is  computed  by 
the  method  of  Ackers  and  White  (1973),  which  uses  a  transport  power  (work 
rate)  approach.  It  has  been  shown  to  provide  superior  results  for  transport 
under  steady-flow  conditions  (White,  Milli,  and  Crabbe  1975)  and  for  combined 
waves  and  currents  (Swart  1976).  Flume  tests  at  the  US  Army  Engineer  Water¬ 
ways  Experiment  Station  have  shown  that  the  concept  is  valid  for  transport  by 
estuarine  currents. 

26.  The  total  load  transport  function  of  Ackers  and  White  is  based  upon 
a  dimensionless  grain  size 


D 

gr 


g(s  -  1) 


1/3 


where 

D  =  sediment  particle  diameter 
s  =  specific  gravity  of  the  sediment 
V  =  kinematic  viscosity  of  the  fluid 
and  a  sediment  mobility  parameter 
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F 

gr 

where 

T  =  total  boundary  shear  stress 

n'  =  a  coefficient  expressing  the  relative  Importance  of  bed-load  and 
suspended-load  transport,  given  in  Equation  A9 

t'  =  boundary  surface  shear  stress 

The  surface  shear  ctress  is  that  part  of  the  total  shear  stress  which  is  due 
to  the  rough  surface  of  the  bed  only,  l.e.,  not  including  that  part  due  to  bed 
forms  and  geometry.  It  therefore  corresponds  to  that  shear  stress  that  the 
flow  would  exert  on  a  plane  bed. 

27.  The  total  sediment  transport  is  expressed  as  an  effective 
concentration 


pgD(s  -  1) 


1/2 
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m 


where  U  is  the  average  flow  speed,  and  for  1  <  <  60 


n  =  1.00  -  0.56  log  D 


gr 


A  =  -2^  +  0.14 


log  C  =  2.86  log  D  -  (log  D  )  -  3.53 

gr  gr 
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(All) 


m  = 


9.66 

D 

gr 


1.34 
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For  D  <60 
gr 


n'  =  0.00 


(A13) 
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A  =  0.17 


CA14) 


C  =  0.025 


(A15) 


m  =  1.5 


(A16) 


28.  Equations  A6-A16  result  in  a  potential  sediment  concentration  G 

P 

This  value  is  the  depth-averaged  concentration  of  sediment  that  will  occur  if 
an  equilibrium  transport  rate  is  reached  with  a  nonlimited  supply  of  sediment. 
The  rate  of  sediment  deposition  (or  erosion)  is  then  computed  as 


where 


C  =  present  sediment  concentration 

t  =  time  constant 
c 

For  deposition,  the  time  constant  is 


I  At 

or 

C  ,h 
d 

V 


and  for  erosion  it  is 


t^  =  larger  of 


At 

or 

C  h 
e 

U 
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where 

At  =  computational  time-step 

=  response  time  coefficient  for  deposition 
Vg  =  sediment  settling  velocity 

*  response  time  coefficient  for  erosion 


AlO 


The  sand  bed  has  a  specified  initial  thickness  which  limits  the  amount  of  ero¬ 
sion  to  that  thickness. 

Cohesive  sediments  transport 

29.  Cohesive  sediments  (usually  clays  and  some  silts)  are  considered  to 
be  depositional  if  the  bed  shear  stress  exerted  by  the  flow  is  less  than  a 
critical  value  .  When  that  value  occurs,  the  deposition  rate  is  given  by 

Krone's  (1962)  equation 


where 

S  =  source  term 

=  fall  velocity  of  a  sediment  particle 
h  =  flow  depth 

C  =  sediment  concentration  in  water  column 

T  =  bed  shear  stress 

=  critical  shear  stress  for  deposition 

C  =  critical  concentration  =  300  mg/£ 
c 


(A20) 


(A21) 


30.  If  the  bed  shear  stress  is  greater  than  the  critical  value  for  par¬ 
ticle  erosion  t  ,  material  is  removed  from  the  bed.  The  source  term  is  then 
e 

computed  by  Ariathural's  (Ariathurai,  MacArthur,  and  Krone  1977)  adaptation  of 
Partheniades '  (1962)  findings: 


S 


for  T  >  T 

e 
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where  P  is  the  erosion  rate  constant,  unless  the  shear  stress  is  also 
greater  than  the  critical  value  for  mass  erosion.  When  this  value  is 
exceeded,  mass  failure  of  a  sediment  layer  occurs  and 


All 
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S  = 


T  P 
L  L 

hAt 


for  T  >  T 

s 


where 

T  =  thickness  of  the  failed  layer 

P  =  density  of  the  failed  layer 

Li 

At  =  time  interval  over  which  failure  occurs 

T  =  bulk  shear  strength  of  the  layer 
s 

31.  The  cohesive  sediment  bed  consists  of  1  to  10  layers,  each  with  a 
distinct  density  and  erosion  resistance.  The  layers  consolidate  with 
overburden  and  time. 

Bed  shear  stress 

32.  Bed  shear  stresses  are  calculated  from  the  flow  speed  according  to 
one  of  four  optional  equations;  the  smooth-wall  log  velocity  profile  or 
Manning  equation  for  flows  alone;  and  a  smooth  bed  or  rippled  bed  equation  for 
combined  currents  and  wind  waves.  Shear  stresses  are  calculated  using  the 
shear  velocity  concept  where 


T 


b 
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where 

T,  =  bed  shear  stress 
b 

=  shear  velocity 


and  the  shear  velocity  is  calculated  by  one  of  four  methods: 
a.  Smooth-wall  log  velocity  profiles 


=  5.75  log 
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which  is  applicable  to  the  lower  15  percent  of  the  boundary 
layer  when 


A12 


>  30 


where  u  is  the  mean  flow  velocity  (resultant  of  u  and  v 
components) 

b.  The  Manning  shear  stress  equation 


CME  (h) 


1/6 


(A26) 


where  CME  is  a  coefficient  of  1  for  SI  (metric)  units  and 
1.486  for  non-SI  units  of  measurement. 

c.  A  Jonsson-type  equation  for  surface  shear  stress  (plane  beds) 
caused  by  waves  and  currents 


u 


* 


(A27) 


where 

f  =  shear  stress  coefficient  for  waves 
w 

u  =  maximum  orbital  velocity  of  waves 
om  ■' 

f  =  shear  stress  coefficient  for  currents 
c 


d.  A  Bijker-type  equation  for  total  shear  stress  caused  by  waves 
and  current 


u 


* 


f  u^ 
w  om 
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Solution  method 

33.  Equation  A5  is  solved  by  the  finite  element  method  using  Galerkin 
weighted  residuals.  Like  RM/.-2V,  which  uses  the  same  general  solution  tech¬ 
nique,  elements  are  quadrilateral  and  may  have  parabolic  sides.  Shape  func¬ 
tions  are  quadratic.  Integration  in  space  is  Gaussian.  Time-stepping  is 


A13 


performed  by  a  Crank-Nicholson  approach  with  a  weighting  factor  (0)  of  0.66. 

A  front-type  solver  similar  to  that  is  RMA-2V  is  used  to  solve  the  simultane¬ 
ous  equations. 
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